
## The following may need to be installed prior to running the code
## install.packages("list") 
## install.packages("multiwayvcov") 
## install.packages('lmtest', dependencies = TRUE)

## RStudio Version 1.1.463 – © 2009-2018 RStudio, Inc. 

# setwd("/Users/Guillem/Downloads/Riambau_Ostwald_Replication_Materials_PSRM/")

library(list)
library(multiwayvcov)
library(lmtest)
sessionInfo()
mydata=load("Data_R_PSLE.rdata")

sink("RIAMBAU_OSTWALD_log_file.txt")

getwd()
set.seed(19791019)   # birthday of one of the authors
options(warn=2)

# checking
mean(Data_R_PSLE$list[Data_R_PSLE$placebo==0])  ## 1.767594
mean(Data_R_PSLE$list[Data_R_PSLE$placebo==1])  ## 1.890555 

#########################
##     CONSTRAINED     ## 
#########################

## PAPER: Without controls 
## [Table 4 S.M., Column 1]
constrained_col1 <- ictreg(list ~ years_school,    data = Data_R_PSLE, treat = "placebo", J=4, method = "ml", overdispersed = FALSE, constrained = TRUE, maxIter = 10000 )
summary(constrained_col1) 
coeftest(constrained_col1)
obs<-lm(list ~    years_school  ,    data = Data_R_PSLE ) # regression run just to confirm number of observations
nobs(obs)  # number of observations confirmed to 1256

## PAPER: With unreported controls (gender, ethnicity, apartment) 
## [Table 4 S.M., Column 2]
constrained_col2 <- ictreg(list ~   female_dummy + female_NA +  apt_3_4_room + apt_5plus_room + Chinese   + years_school,
                           data = Data_R_PSLE, treat = "placebo", J=4, method = "ml", 
                           overdispersed = FALSE, constrained = TRUE, maxIter = 10000 )
summary(constrained_col2)
coeftest(constrained_col2)
obs<-lm(list ~ female_dummy + female_NA +  apt_3_4_room + apt_5plus_room + Chinese   + years_school,    data = Data_R_PSLE ) # regression run just to confirm number of observations
nobs(obs)  # number of observations confirmed to 1254 


## PAPER: With reported controls (age, income)
## [Table 4 S.M., Column 3]
constrained_col3 <- ictreg(list ~   years_school + age_61plus + age_NA + hhd_income_in_thousands   ,    data = Data_R_PSLE, treat = "placebo", J=4, method = "ml", 
                                                overdispersed = FALSE, constrained = TRUE, maxIter = 10000 )
summary(constrained_col3)
coeftest(constrained_col3)
obs<-lm(list ~ years_school + age_61plus + age_NA + hhd_income_in_thousands,    data = Data_R_PSLE ) # regression run just to confirm number of observations
nobs(obs)  # number of observations confirmed to 1158 

## PAPER: With reported and unreported controls (age, income, gender, ethnicity, apartment)
## [Table 4 S.M., Column 4]





  # Error in solve.default(-MLEfit$hessian) : 
  # system is computationally singular: reciprocal condition number = 3.55341e-2

  ## Using "car" as a proxy for income instead of apartment size  (else the code will not converge)

  constrained_col4 <- ictreg(list ~ years_school + age_61plus + age_NA +  hhd_income_in_thousands
                           + female_dummy + female_NA + car + Chinese,
                           data = Data_R_PSLE, treat = "placebo", J=4, method = "ml", 
                           overdispersed = FALSE, constrained = TRUE, maxIter = 10000 )
  summary(constrained_col4)
  coeftest(constrained_col4)
  obs<-lm(list ~ years_school + age_61plus + age_NA +  hhd_income_in_thousands +
          female_dummy + female_NA + car + Chinese,    data = Data_R_PSLE ) 
  # regression run just to confirm number of observations
  nobs(obs)  # number of observations confirmed to 1156
 
  
#########################
##    UNCONSTRAINED    ## 
#########################

## PAPER: Without covariates, UNCONSTRAINED
## [Table 4 S.M., Column 5]
            # Changing warning options (else converging issues)
            options(warn=1)
            unconstrained_col5 <- ictreg(list ~ years_school, data = Data_R_PSLE,
                                  treat = "placebo", J=4, method = "ml",
                                  overdispersed = FALSE, constrained = FALSE, maxIter = 10000)
            summary(unconstrained_col5)
            coeftest(unconstrained_col5)
            # Changing warning optionsback to "2"
            options(warn=2)
 
## PAPER: With unreported controls (gender, ethnicity, apartment) 
## [Table 4 S.M., Column 6]
  
## Using "car" as a proxy for income instead of apartment size (else the code will not converge)
  unconstrained_col6 <- ictreg(list ~  years_school +  female_dummy  + female_NA + Chinese + car,
                             data = Data_R_PSLE, treat = "placebo", J=4, method = "ml", 
                             overdispersed = FALSE, constrained = FALSE, maxIter = 10000)
  summary(unconstrained_col6)
  coeftest(unconstrained_col6)
 
## PAPER: With reported controls (age, income)
## [Table 4 S.M., Column 7]
## Yields intercept much bigger than 5 (8.15!)
unconstrained_col7 <- ictreg(list ~ years_school + age_61plus + age_NA + hhd_income_in_thousands   , data =Data_R_PSLE, treat = "placebo", J=4, method = "ml", 
                                                overdispersed = FALSE, constrained = FALSE, maxIter = 10000)
summary(unconstrained_col7)
coeftest(unconstrained_col7)
 
## PAPER: With reported and unreported controls (age, income, gender, ethnicity, apartment)
## [Table 4 S.M., Column 8]
## Using "car" as a proxy for income instead of apartment size (else the code will not converge)
  unconstrained_col8 <- ictreg(list ~ years_school + age_61plus + age_NA + hhd_income_in_thousands
                             + female_dummy + female_NA  +  Chinese  + car,
                             data = Data_R_PSLE, treat = "placebo", J=4, method = "ml", 
                             overdispersed = FALSE, constrained = FALSE, maxIter = 10000)
  summary(unconstrained_col8)
  coeftest(unconstrained_col8)
 
  
 